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|  ABSTRACT 
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The  specific  impetus  for  this  work  was  a  conceptual  design  of  a 
submarine  using  the  toroid  as  the  pressure  hull.  This  work  is  a 
continuation  of  the  work  started  by  Bowen  (1987).  As  such,  it  is  hoped 
that  a  better  understanding  of  the  behavior  of  a  toroid  under 
hydrostatic  pressure  can  be  realized. 

This  work  began  with  a  review  of  efforts  to  solve  complete  toroidal 
structures.  A  specific  toroid  was  then  modeled  in  the  B0S0R4  computer 
program  to  obtain  displacements  of  the  meridian  under  hydrostatic 
pressure.  Functions  were  then  derived  that  described  the  general  form 
of  these  displacements.  Using  these  functions  as  the  assumed  solution 
for  the  energy  method  an  energy  balance  was  made  and  a  program  was 
written  to  solve  for  the  displacements  of  a  generic  toroid  under 

hydrostatic  loading.  K  ,  \  f, 
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CHAPTER  1 
INTRODUCTION 


The  purpose  of  this  work  was  to  examine  the  response  of  a  thin 
circular  toroidal  shell  under  external  hydrostatic  pressure.  The 
approach  used  was  to  model  a  specific  toroid  in  the  B0S0R4  computer 
program  (Bushnell,  1977)  and  then  use  the  output  results  for 
displacements  to  generate  functions  that  modeled  these  displacements. 
Using  the  functions  obtained  that  model  the  displacements  of  the  toroid 
under  hydrostatic  pressure,  an  analytical  analysis  was  performed  with 
these  functions  used  as  the  assumed  solution  for  the  energy  method.  The 
energy  method  generated  the  equations  necessary  to  solve  for  the 
displacements  of  a  generic  circular  toroid.  A  program  was  then  written 
to  solve  the  equations  generated  in  the  energy  method  so  that 
displacements  for  any  specific  toroid  under  hydrostatic  pressure  could 
be  obtained. 

The  analysis  conducted  only  considers  displacements  in  the  linear 
elastic  range.  For  the  shell  being  analyzed  the  following  assumptions 
were  made : 

1.  The  material  of  the  shell  is  isotropic  and  homogeneous. 

2.  The  thickness  of  the  shell  is  constant. 

3.  The  thickness  of  the  shell  is  small  compared  to  the  radii  of 

curvature . 

b.  The  displacements  are  symmetric  about  the  X-Y  plane  (see 

Figure  1.1). 

5.  Thin  shell  theory  holds:  normals  to  the  undeformed  surface 
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remain  normal. 


For  a  symmetrically  loaded  shell,  with  the  above  assumption  of 
symmetric  deformation,  a  small  displacement  of  any  point  can  be  resolved 
into  two  components:  v  in  the  direction  of  the  tangent  to  the  meridian 
and  u  in  the  direction  of  the  normal  to  the  middle  surface  (Timoshenko 
and  Woinowsky-Krieger ,  1959).  Therefore,  the  only  displacements 
discussed  throughout  the  text  will  be  u  and  v. 

The  toroid's  geometry  is  straight  forward  but  in  order  to  be 
consistent  throughout  the  text  the  following  definitions  will  be  made 
and  are  illustrated  in  Figure  1.1: 


DEFINITIONS : 

^  R  :  Major  radius  of  toroid.  Radius  of  rotation  about  the  Z  axis 

of  the  circle  of  radius  r  to  form  the  toroid. 

r  :  Minor  radius  of  toroid.  Radius  of  the  circle  which  is  rotated 
about  the  Z  axis . 

I 

h  :  Thickness  of  the  shell  (assumed  to  be  constant) . 

8  Angle  of  rotation  measured  counter  clockwise  from  the  X-Y 
plane  at  a  distance  R  +  r  from  the  origin. 

i 

<f>  :  Angle  of  rotation  about  the  Z  axis  measured  counter  clockwise 
from  the  positive  X  axis . 


i 


t 
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Figure  1.1 


Description  of  Geometry 
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CHAPTER  2 
COMPUTER  MODELING 

A  computer  solution  for  the  displacements  of  the  toroid  under 
external  hydrostatic  pressure  was  sought  to  determine  the  shape  of  these 
displacements  for  a  specific  toroid.  The  computer  program  used  to  model 
the  toroid  was  B0S0R4  (Bushnell,  1977).  The  B0S0R4  computer  program  is 
a  hybrid  finite  difference  -  finite  element  program  used  to  analyze 
complex  shells  of  revolution.  This  program  was  developed  by  David 
Bushnell  at  Lockheed  Missiles  &  Space  Co.,  Palo  Alto,  California. 

Although  specific  reference  was  not  found  that  indicated  that 
B0S0R4  could  adequately  model  a  tor>id,  it  was  believed  that  since 
B0S0R5 ,  a  similar  program  that  considers  elastic -plastic  material 
behavior,  had  been  used  to  model  a  toroid  that  B0S0R4  would  also 
adequately  model  this  structure  (Bushnell,  1985). 

The  specific  toroid  that  was  analyzed  with  B0S0R4  was  the  toroid 
specified  as  the  default  values  by  Bowen  (1987) .  The  toroid  geometry 
used  was  as  follows: 

R  -  8  inches 
r  -  2  inches 
h  -  0.02  inches 

Young's  Modulus,  E  -  30  x  106  psi 
Poisson's  Ratio,  v  -  0.3 

Modeling  the  toroid  in  B0S0R4  was  accomplished  by  taking  advantage 
of  the  symmetry  of  loading  of  the  structure  about  the  X-Y  plane. 
Therefore,  only  half  of  the  structure  was  discretized  in  the  program. 


Because  only  half  of  the  structure  was  modeled  in  the  program,  boundary 
conditions  had  to  be  established  at  the  symmetry  points.  The  boundary 
conditions  imposed  on  the  structure  are  as  follows: 


@ 


@ 


8-0  radians 

Tangential  Displacement,  (v)  - 
Rotation  in  the  9  direction, 

9  -  II  radians 

Tangential  Displacement,  (v)  - 
Rotation  in  the  9  direction, 


0.0 
-  0.0 


0.0 
-  0.0 


The  normal  displacement  (u)  is  defined  as  positive  in  the  direction 
of  the  outward  pointing  normal.  The  tangential  displacement  is  defined 
as  positive  in  the  direction  of  increasing  angle  9  (See  Figure  2.1). 

Having  the  toroid  modeled  in  B0S0R4,  the  program  was  run  to 
determine  the  pressure  required  to  buckle  the  structure .  The  thought 
being  that  if  the  buckling  pressure  predicted  by  B0S0R4  was  consistent 
with  the  buckling  pressures  predicted  by  theory  then  the  pre -buckling 
displacements  generated  by  B0S0R4  could  be  used  to  predict  the  general 
form  of  the  displacements  of  a  generic  toroid. 

The  buckling  pressure  obtained  from  B0S0R4  for  the  toroid  modeled 
differed  by  less  than  10%  from  the  predicted  buckling  pressure  as 
presented  by  Sobel  and  Flugge  (1967).  Since  the  buckling  pressure 
obtained  was  consistent  with  the  predicted  pressure  and  the  predicted 
buckling  pressure  was  substantiated  experimentally  by  Almroth,  Sobel  and 
Hunter  (1969),  it  was  concluded  that  the  results  from  B0S0R4  for  the 
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Figure  2.1 


B0S0R4  Modeling 


pre -buckling  displacements  could  be  used  to  predict  the  general  form  of 
the  displacements  of  a  generic  toroid  loaded  under  uniform  hydrostatic 
pressure . 

In  any  finite  difference  or  finite  element  program  an  investigation 
into  the  number  and  spacing  of  the  nodes  used  to  describe  the  structure 
must  be  undertaken  to  ensure  that  the  structure  is  adequately  described. 

For  this  investigation,  equal  spacing  of  the  nodes  was  used  throughout. 

To  investigate  the  effect  of  node  distribution,  the  number  of  node:, 
were  doubled  and  thus  the  spacing  was  reduced  by  a  factor  of  two.  The 
results  for  both  the  magnitude  of  the  buckling  pressure  and  the  shape 
and  magnitude  of  the  displacements  did  not  change  with  the  increase  in 
the  number  of  nodes.  It  was  therefore  concluded  that  the  structure  was 
adequately  described  and  that  the  results  obtained  were  valid. 

Figure  2.2  is  a  graphical  presentation  of  the  normal  displacement 
(u)  and  figure  2.3  is  a  graphical  presentation  of  the  tangential 
displacement  (v)  obtained  from  the  B0S0R4  output  for  the  toroid  modeled. 
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CHAPTER  3 


CURVE  FITTING 

The  next  step  in  this  study  was  to  generate  a  function  that  would 
analytically  describe  the  normal  and  tangential  displacements  of  the 
toroid  under  hydrostatic  loading. 

First  a  check  was  made,  using  the  B0S0R4  computer  program,  to 
ensure  the  displacements  were  linearly  dependent  on  the  applied 
pressure,  as  should  be  the  case  since  this  study  was  only  considering 
behavior  of  the  toroid  in  the  linear  elastic  region.  As  can  be  seen 
from  figure  3.1  the  normal  displacement  is  indeed  linearly  dependent  on 
the  applied  pressure.  Although  not  included  here,  the  tangential 
displacement  is  also  linearly  dependent  on  the  applied  pressure. 

To  generate  a  function  that  would  adequately  describe  the 
displacements,  a  fora  of  the  displacements  must  be  assumed  and  then  a 
curve  fitting  technique  must  be  used  to  fit  the  data  to  this  assumed 
form. 

Since  the  displacements  are  considered  to  be  symmetric  about  the 
X-Y  plane  and  considering  the  way  in  which  the  positive  direction  for 
the  normal  displacement  (u)  and  the  tangential  displacement  (v)  were 
defined  in  Chapter  2,  the  following  statements  can  be  made: 

For  the  normal  displacement  (u) , 

u(F)  -  u(-f)  (3.1) 

which  is  the  definition  of  an  odd  function. 

For  the  tangential  displacement  (v) , 

v($)  -  -v(-F)  (3.2) 
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Figure  3.1 
Increment*!  Loading 


Formal  Displacement 


which  is  the  definition  of  an  even  function. 

The  functions  assumed  to  describe  the  displacements  were  Fourier 
Series  for  both  the  normal  (u)  and  tangential  (v)  displacements.  Since 
both  the  normal  and  tangential  displacements  have  a  period  of  2irr  (the 
arc  length  of  the  circular  cross  section)  and  noting  the  above 
statements  about  even  and  odd  functions,  the  Fourier  Series  for  the 
displacements  can  be  written  as  follows: 

u(0)  -  ~  +  X  an  cos  (3.3)  , 

n-l  ^  * 

V(t)  -  £  bn  sin  (3.4). 

For  the  displacements  to  be  symmetric  about  the  X-Y  plane,  or  a 
half  period  of  e  function,  x  and  L  can  be  defined  as  follows: 

x  -  r  9 
L  -  irr 

From  these  definitions,  u  and  v  can  be  rewritten  as: 


U($)  -  cos  (n0) 

n-l 

00 

v(0)  -  JT  bn  sin  (nf) 
n-l 


(3.5), 


I 


(3.6). 


While  the  Fourier  series  described  above  will  give  an  exact 
representation  of  the  displacements,  the  goal  of  a  curve  fitting  routine 
should  be  to  determine  how  many  terms  of  the  series  must  be  used  to 


I 
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adequately  describe  the  function.  To  determine  the  number  of  terms 
required  to  describe  the  displacements  a  least  squares  curve  fitting 
program  was  written  to  be  used  in  conjunction  with  the  output  data  from 
B0S0R4.  The  programs  were  adapted  from  a  curve  fitting  routine 
presented  by  James,  Smith  and  Wolford  (1977)  and  are  included  in 
Appendix  A. 

Using  these  programs  to  determine  the  coefficients  of  the  Fourier 
Series  it  was  found  that  five  terms  were  inadequate  to  describe  the 
displacements,  as  can  be  seen  graphically  in  Figure  3.2,  and  it  was 
determined  that  ten  terms  were  required  to  adequately  describe  the 
normal  displacements,  (See  Figure  3.3),  and  nine  terms  were  required  to 
describe  the  tangential  displacements.  Although  only  the  normal 
displacement  (u)  is  presented  in  Figures  3.2  and  3.3  similar  results 
were  obtained  for  the  tangential  displacement  (v) . 

With  the  number  of  terms  determined  that  are  required  in  the 

Fourier  series  u  and  v  can  be  written  in  final  form  as: 

9 

u(0)  -  +  X  an  cos  (n 9)  (3.7), 

n  •  1 
9 

v(0)  “  £  bn  sin  (n0)  (3.8). 

n«l 

These  definitions  of  u  and  v  will  be  used  in  the  next  section  as  rhe 
assumed  forms  of  the  displacements  of  the  toroid. 
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CHAPTER  4 


ENERGY  METHOD 

4.1  Introduction: 

To  determine  the  displacements  of  a  structure  when  it  is  subjected 
to  an  external  load,  the  external  load  must  be  balanced  by  the  internal 
forces  of  the  structure  caused  by  the  external  load. 

For  this  study,  the  method  of  minimizing  the  total  potential 
energy,  as  described  by  Shames  and  Dym  (1985),  was  used.  The  energy 
method  was  cho  en  in  hopes  of  avoiding  the  singularities  in  the  solution 
at  8  -  ±  -  .(Reissner,  1963  and  Senjanovic  ,  1972).  The  total  potential 
energy  is  defined  as: 

n  -  u  -  w  (4.1.1) 

where  U  and  W  are  defined  as; 

U  :  Internal  energy  of  the  structure, 

W  :  potential  energy  of  the  applied  loads. 

To  determine  the  minimum  potential  energy,  and  thereby  determine  the 
displacements  caused  by  the  external  load,  the  first  variation  of  the 
total  potential  energy  is  set  to  zero. 

6U,n  -  0  (4.1.2) 


The  first  variation  of  the  total  potential  energy  can  be  defined  as  the 
following  partial  differential  equation: 


an 

3x 


-  0 


(4.1.3) 


i 

where  is  any  general  parameter  used  to  describe  U  and  W. 

Since  the  first  variation  of  the  total  potential  energy  is  set  to  zero, 


I 


I 


I 


« 
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Che  equation  Co  be  solved  to  determine  the  displacements  of  the 


structure  can  be  defined  as; 


5W  _  5U 
dx  3x 


(4.1.4) . 


In  order  to  solve  the  above  equation  for  the  displacements  caused  by  the 
applied  load,  U  and  W  must  be  defined  as  functions  of  the  disj.  .acements . 
As  described  by  Bushnell  (1984),  the  strain  energy  of  the  structure 


can  be  defined  as  follows: 


where ; 


u  -  u  +  u 

s  c 


(4.1.5) 


Us  strain  energy  due  to  elongation  and  changes  in 
curvature , 

U_  :  strain  energy  due  to  constraints  or  boundary  conditions. 
For  the  toroid,  the  symmetry  about  the  X-Y  Plane  will  be  taken 
advantage  of  and  only  half  of  the  structure  will  be  analyzed.  The 
following  additional  geometric  definitions  are  required  for  the  analysis 
(See  Figure  4.1). 


DEFINITIONS: 


Ratio  of  the  major  and  the  minor  radius  of  the  toroid,  a  -  - 
Distance  from  the  Z  axis  to  the  meridian  of  the  toroid, 
r^  -  R  +  r  cos(fl)  -  r[  a  +  cos  (0)) 

Radius  of  curvature  in  the  9  direction. 


t  _  R  +  r  cos ( 9 ) 
ra  “  cos(0)  ”  cos(0) 
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4.2  Strain  Energy: 

The  strain  energy  due  to  elongation  and  changes  in  curvature  can  be 
defined  as: 

U  -  —  T  a  £  dv  ,  i.j  -  9  ,  4  (4.2.1). 

S  2  Jv  tj  ij 

The  strain  energy  due  to  constraints,  U  ,  will  be  addressed  later. 

Because  this  study  is  only  considering  hydrostatic  loading  of  the 
toroid,  the  loading  will  always  be  normal  to  the  surface  of  the  toroid. 
Therefore,  there  will  not  be  any  shear  stresses  or  shear  resultants  from 
the  loading  imposed. 

a  —  f  »  0  ,  if  i  »<  j. 

U  tj 

The  8  and  4  directions  are  thus  principle  directions  and  the  strains  can 
be  represented  by  reference  surface  strains  and  changes  in  curvature. 
The  strain  energy  can  be  defined  in  terms  of  forces  and  moments  which 
are  defined  as  an  integral  of  the  stresses  through  the  thickness  of  the 
shell.  The  strain  energy  can  then  be  represented  by  a  surface  integral 
as  follows: 

U-if  f  N  e  +  M  *  )ds  (4.2.2). 

s  2  JI  *■  U  ij  tj  ij  J 

For  the  case  of  hydrostatic  loading  the  strain  energy  becomes: 

u.  -  jj.tcv, +  v*)  *  (V» 4  V*  )]ds  “'2  3) 

where  the  forces  and  moments  are  defined  as  follows  (Timoshenko  and 
Uoinowsky  -  Krieger,  1959): 

N*  ■  jt  °#  (  1  '  I  ]  dz 

- 1.  %  [ 1  '  fj  d2  (‘-2'5>' 
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*9  ~  Jt  2  C  1  *  f  )  dz  (4.2.6) , 

-  I,  %  2  (  1  -  ; J  dz  <4-2-7>- 

Since  the  assumption  has  been  made  that  normals  to  the  undeformed 
middle  surface  remain  normal,  the  strains  as  a  function  of  thickness  can 
be  expressed  in  terms  of  reference  surface  strains  and  changes  in 
curvature  (Timoshenko  and  Woinovsky-Krieger , 1959) . 

e$  -  «9  ■  z  k6  (4.2.8) 


-  t 


-  Z  K 


(4.2.9) 


t 

The  terms,  e  and  «  ,  represent  the  total  strain  in  the  9 

u  © 

t  t 

directions  and  k  and  k,  represent  the  changes  in  curvature 
reference  surface. 

From  Hook's  Law  the  stresses  can  be  defined  as  follows: 


and  4> 
of  the 


%  " 


"  -1—  f  ^  1 

(1-y2)  9t 

h:  ( '*  * v<* ) 

-v  )  ^  t  tr 


(1*0 

substituting  the  expressions  for  t.  and  <  , 

t  t 

’•  -  ^77  ( e» 4  v'*  ■<  "> +  “% )  ] 
°* '  7^7")  ( '* +  'zC  +  *"•* 5 ) 


(4.2.10) 


(4.2.11) 


(4.2.12), 


(4.2.13) . 


Taking  the  middle  surface  as  the  reference  surface  and  neglecting 
z  z 

the  terms,  -  and  -  ,  as  small  as  compared  to  unity,  integration  over 
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I 

I 

I 

I 

I 

I 


> 


u +  ue* 


the  thickness  of  the  shell  can  be  represented  by: 

r+h/2  E  f  1 

Nf  :J,  (iv)  I  '•  *  "*  'zC  ■* +  ]  J 

N«  ’  7^7  ( *•  *  “'♦ ) 

r*h/2  e  r  1 

♦  J*  (1-uV  *  •  L  *  ,JJ 

N,  --3l-  r 

*  (l-v!)  l 

r+h/2  c  f 

-  J  7—7:  l‘(  * 

-h/2  (1-U  ) 

r+h/2  e  f  1 

M,  -  -  «,  +  vt.  -zfx,  +  VK.  )  z  dz 

*  -J,Z  <lV)  t.  *  >  1  *  >’) 

;h*  f  1 

,-y2)l  *  9  J 


Kg  +  II  Z  dz 


) 

3) 


K.  +  UK. 

9  4> 


) 


-  Eh 


12(1- 


(4.2.14)  , 

(4.2.15)  , 

(4.2.16)  , 

(4.2.17)  , 

(4.2.18)  , 

(4.2.19)  , 

(4.2.20) , 

(4.2.21)  . 


Combining  the  expressions  for  the  forces  and  the  moments  and 
substituting  these  back  into  the  expression  for  the  strain  energy,  the 
strain  energy  becomes: 


U 


1 

2 


Eh 

(lV) 


+  2uV* 


Eh" 


12(1-0 


C  K)  +  2vKeKt  +  % 


3) 


ds 


(4.2.22) . 


4.2.1  Strains: 

The  next  step  is  to  define  the  surface  strains  in  terms  of  the 
displacements  u  and  v.  This  will  be  accomplished  by  following  the 
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presentation  by  Timoshenko  and  Uoinowsky-Krieger , (1959) . 


Considering  an  element  AB  of  the  meridian,  (See  Figure  4.2),  the 
strain  in  the  9  direction  can  be  defined  as  follows: 

«  due  to  v , 


€  “ 


<9v 

V  +  -  d»  -  V  J  ^ 

r  d 8  r  30 


(4. 2. 1.1) 


<  due  to  u  can  be  represented  by  the  average  change  in  length  of 


the  element  divided  by  the  original  length, 

1  5u 

2  88 

(  —  - 

l 

l  du 


(4.2. 1.2) . 


r  d  8 

The  component,  jj  dff  ,  will  be  neglected  as  a  small  quantity  of  higher 
order.  The  total  strain  in  the  6  direction  can  then  be  represented  as: 


'8 


If  5v  1 

l“*n 


(4. 2. 1.3) . 


The  strain  in  the  4  direction  can  be  defined  as  the  change  in  r„ 
due  to  u  and  v  divided  by  the  original  length  of  the  element  (See  Figure 
4.3)  . 

A  r^  -  [  u  cos (8)  -  v  sin(0)  )  d^  ,  and 
[  u  cos (0)  -  v  sin(0)  )  <14 


«  ,  - 


r  d 4 

T 


‘  r  (7  Kos(«))(u  cos(<?)  '  v  sin((?) 


(4. 2.1. 4) , 
(4. 2.1. 5) . 


4.2.2  Change  in  Curvature : 

*  The  change  in  curvature  of  the  middle  surface  will  now  be  defined. 

Again  considering  the  element  AB  of  the  meridian,  (See  Figure  4.2), 

the  rotation  in  the  9  direction  can  be  defined  as  follows: 

I 


1 


p  due  to  v, 


p  -  -  dff 
r 


(4. 2. 2.1) 


P  due  to  u, 


&  - 


3U  AD 

39  d6 


-  1 


3u  At) 
39  d6 


The  total  rotation  in  the  9  direction  is; 


**  -  \[ v  •  It  )  d} 


(4. 2. 2. 2) . 


(4.2.2. 3) 


Th.3  change  in  curvature  in  the  8  direction  can  be  represented  by  the 
change  in  rotation  in  the  9  direction  divided  by  the  undeformed  arc 
length  in  the  8  direction. 


1  ffl 

Ke  m  r  ae 


l  f  3v  a2 u  ) 


(4. 2. 2. 4) 


Because  the  assumption  has  been  made  that  the  deformations  will  be 
symmetric,  the  rotation  in  the  d  direction  will  be  opposite  to  the 
rotation  in  the  9  direction  in  magnitude  and  rotated  into  the  d  surface 
by  sin  (9)  . 

^  sin  (9)  -  |(  fj  -  v  )  sin  ( 8 )  d d  (4. 2. 2. 5) 

Although  the  rotation  in  the  d  direction  is  constant  in  d  ,  it  does  vary 
in  the  9  direction.  Therefore,  the  change  in  curvature  due  to  loading 
can  be  represented  by  the  rotation  in  the  d  direction  divided  by  the 


undeformed  arc  length  in  the  d  direction. 

0, 


K  .  - 


r  (a  +  cos (0) )  dd 


(4. 2. 2. 6) 
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■ 


and  substituting  the  expression  for  fi  , 

sin(fl)  f  3u  1 

ka  "  2  I  38  '  V  I 

9  r  (a  +  cos (9))  ^  J 


(4. 2. 2. 7)  . 


4.2.3  Constraint  Energy : 

A  complete  toroid  has  no  constraint  or  boundary  conditions  when 
loaded  only  by  hydrostatic  pressure.  However,  since  only  half  of  the 
toroid  is  being  analyzed  here,  boundary  conditions  must  be  imposed  at 
the  symmetry  points,  9-0  and  9  -  *.  The  boundary  conditions  that  are 
required  at  these  two  points  are: 

The  vertical  displacement  is  zero,  v  -  0, 

and  the  rotation  in  the  9  direction  is  zero,  ^  -  0. 

0  v 

As  Bushnell  (1984)  points  out  the  constraint  energy  can  be  accounted  for 
with  the  introduction  of  Lagrange  multipliers  and  can  be  represented  as 
follows : 


u  -  +  *  v<°>  +  *  i?<*>  +  *  v(*> 

e  Quo  0  Ov  Xua  v  XV 


(4.2.3.1)  . 


In  this  formulation  of  the  constraint  energy  the  Lagrange  multipliers 
are  unknowns  and  are  solved  for  along  with  the  displacements  when  the 
total  potential  energy  is  minimized.  The  actual  form  of  the  constraint 
energy  will  be  discussed  later. 


4.3  Potential  Energy  of  Applied  Loads: 

The  potential  energy  of  the  applied  loads  is  simply  the  work  done 
by  these  loads,  which  can  be  represented  by  the  applied  load  multiplied 
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by  the  deflection  in  the  direction  of  the  applied  load.  For  the  case  of 
hydrostatic  loading,  the  applied  load  is  always  normal  to  the  surface 
and  the  displacement  in  the  direction  of  this  load  is  -u  for  this  study. 

Therefore,  the  potential  energy  of  the  applied  hydrostatic  load  can  be 
represented  as : 

W  -  -  P  u  ds  (4.3.1) 

where  P  is  the  external  hydrostatic  pressure  applied  and  is 
positive  in  the  direction  of  -u. 


4.4  Minimization  of  the  Total  Potential  Energy: 

From  section  4.1  it  was  shown  that  to  minimize  the  total  potential 
energy  the  following  equation  needed  to  be  solved; 


the 


g-o 

l 

and  it  was  shown  that  the  above  equation  could  be  represented  by 
following  equation: 


aw  _  au 

3x  3x 


(4.4.2)  . 


l  l 

To  present  the  form  of  the  above  equations  the  following  simplifications 
in  notation  will  be  used: 


,  du 

U  “  36  ’  U 


a2u 

,  ,  _  ,  3v 

‘  a*2  ’  "  ' 


c  - 


Eh 


(1-v  ) 


D  - 


Eh' 


12  ( 1  -  u  ) 
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4.4.1  Internal  Energy : 


Since  ic  has  been  shown  that  none  of  the  terms  in  the  internal 
energy  depend  on  4  the  surface  integral  can  be  rewritten  and  integrated 
with  respect  to  6. 

jk  ,zx , 


— — — r  +  2vk.k,  +  e  +  cos(0)]d^  rd 9 

12(lV)  9  9  * 


(4.4. 1.1) 


Which  when  integrated  with  respect  to  ^  becomes, 


U  -  t  r 
8 


f  pV  C 

i  '•(l-u  ) 


ee  +  2uV*  +  0  ' 


12 


^7  ("J  +  2uV*  +  %)](  a  +  C0S^))dJ 


(4. 4. 1.2) . 


Using  the  above  simplifications  in  notation  and  substituting  the 
expressions  presented  in  section  4.2.1  and  4.2.2  the  strain  energy 
becomes : 


U 


x  r 


f(H^ 


+  2uv' 


+  v 


,2 


+ 


(a  +  cos (8) )  '• 

+  _ l _  2 

(a  +  cos (9)) 

v2sin 2 (9)  )  j  - 

2  v  sin(fl)  r 
(a  +  cos(  9 ) )  *■ 


(u2+uv'  )cos(0 )  -  (uv  +  w')sin(0)) 
[  u2cos2(0)  -  2uv  cos ( 9 ) sin( 9 )+ 

^  (v’2-  2v'u”+  u"2  )  + 
v'u'-  w '  -u'u' '  +vu ' ' ]  + 
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(«  2“'v  *  v1]]]  C.«o.<.))« 


(4.4. 1.3) . 


The  above  equation  represents  the  strain  energy,  U  .  However,  for 

minimization  of  the  total  potential  energy  the  expression  that  is  needed 

.  3U 
is  -7—  . 

3x 

i 

au 


ar  -  *  I  ( c  ( C  ^  *  *•£■) 

i  ii  t 


2u  r/r.  <9u  dv‘  rdu  x  .  dv  du  £v'  ,  £v  «v 

U~~o.T«n((2“H  3x  )cos(().(uH  ♦*£  *v  _  )sin(tf )) 


'(«•-  c.«(»))aP^;  os!«)-2(u£  )co.(«),ln(»)*2v|i  .ln2<»>]] 

1  11  1  1 

■  [  C2v'h  •  2<v'l;';  “"h  ’-  2u"Ij")  * 

'•lit  l 

2  i/  sin($)  r  ,3u'  ,3v'  3v'  ,3v  ,au''  ,  .du'  3u'  ’  ,,3v-> 

(a  +  cos(a))1-  ax  ax  ax  ax  ax  ax  ax  ax  J 

liii  i  i  i  i 

u~‘1cos(7))iC2u'l;  -2<u'M  +vh )+2vfx 3]) [*«<-•<*»"  <‘■'■•1 

l  ii  l 

The  above  equation  represents  the  first  term  on  the  right  hand  side  of 
the  following  equation  that  will  be  needed  to  minimize  the  total 
potential  energy  and  thereby  solve  for  the  displacements, 

aw 


au  au 

t  c 

ax  ax  ax 

i  i  i 


(4. 4. 1.5) . 


4.4.2  External  Energy : 

Similar  to  the  internal  energy,  the  external  energy  does  not  depend 
on  $  and  the  surface  integral  can  be  rewritten  and  integrated  with 
respect  to  4>  ■ 


“  -  •  rn  P  u  j  r[  a  +  cos(S)  ]  d^  r  d 9 


(4. 4. 2.1) 
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Which  when  integrated  with  respect  to  becomes, 


r* 

W  -  -2P*r2  |  u  (  a  +  cos(0)  )  d8 


(4. 4. 2. 2) 


For  minimization  of  the  total  potential  energy  the  expression  that  is 

.  .  .  aw 

needed  is  ^ 

i 

|^  -  -2Pirr2  J  (  a  +  cos(*)  ]  d6  (4. 4. 2. 3) 

i  u  1 

This  completes  the  expressions  needed  for  the  minimization  of  the  total 
potential  energy,  with  the  exception  of  the  constraint  energy,  U 


4.5  Assumed  Functions: 

In  Chapter  3  it  was  shown  that  the  displacements  u  and  v  could  be 
represented  by  a  Fourier  series.  These  functions  of  9  will  be  used  as 
the  assumed  form  of  the  displacements  and  the  coefficients  of  the 
Fourier  series  will  be  solved  for  from  the  equations  presented  in  the 
minimization  of  the  total  potential  energy.  The  assumed  functions  for  u 
and  v  are  different  from  those  presented  in  Chapter  3  in  that  the 
constant  term  in  the  expression  for  u  has  been  included  inside  the 
summation. 

9 

u(5)  -  \  an  cos  (nfl)  (4.5.1), 

n»0 

9 

v(0)  -  £  bn  sin  (n 6)  (4.5.2). 

n«0 

The  above  expressions  for  u  and  v  involve  19  unknown  coefficients  that 
will  have  to  be  solved  for  to  describe  the  displacements . 

Using  these  expressions  for  u  and  v  the  constraint  energy,  , 


discussed  in  section  4.2.3  can  now  be  addressed.  The  constraints  needed 
were  the  boundary  conditions  at  the  symmetry  points,  9-0  and  9  -  n. 
These  boundary  conditions  are: 

v  -  0.  Since  sin(ntf)  -  0,  for  9-0  and  sin(n»r)  -  0,  for  all  n, 
these  boundary  conditions  are  satisfied  by  the  assumed  form  of  v 
and  additional  constraints  are  not  needed. 

9 

|^-0.  Since  -  -  Y  n  as  sin  (n 9)  ,  then  ^  •  0  at  }  *  0  and 
off  off  off 

n  *0 

9  -  n  for  the  same  reasons  as  stated  above  for  v. 

Since  no  additional  constraints  are  needed  to  meet  the  boundary 
conditions,  Lagrange  multipliers  are  not  needed  and  the  constraint 
energy,  U  -  0. 


4.6  Solution: 

With  the  assumed  forms  of  the  displacements  defined,  the  problem  is 
now  restricted  to  solving  for  the  19  unknown  coefficients  of  the  Fourier 
series.  Using  the  following  substitutions  the  equations  can  be  set  up 


to  solve  for  these  coefficients: 

9 

u-  £ancos(n0)  (4.6.1), 

n  *0 

9 

u'  -  -  ^  n  an  sin  (n 9)  (4.6.2), 

n  *0 
9 

u' '  -  -  £  n2  an  cos  (n0)  (4.6.3), 

n  *0 
9 

v  -  bn  sin  (n£)  (4.6.4), 

n  «0 
9 

v'  -  £  n  bn  cos  (n9)  (4.6.5). 

n  *0 
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Substituting  the  above  expressions  into  the  following  equation  will  then 
represent  the  19  equations  needed  to  solve  for  the  19  unknown 
coeff icients : 


dx  dx 

i  i 


Which  can  be  represented  in  matrix  notation  as  follows; 


C  D  ]  -  [  F  ]  [  x  ] 


(4.6.6) . 


(4.6.7) 


where  the  above  terms  are  defined  as: 


D  :  a  19  x  1  column  consisting  of  the  terms,  —  , 

ox 

l 

3U 

F  ;  a  19  x  19  matrix  consisting  of  the  terms,  — *  , 

i 

x^  :  a  19  x  1  column  consisting  of  the  terms,  an  and  bn.  Where  x^ 
through  x  are  the  10  an's  and  x  through  x^g  are  the  9 
bn's. 

To  solve  for  the  x  's,  the  coefficients  of  the  Fourier  series,  the 

i 

following  matrix  operation  needs  to  be  performed: 


[\i  -  c*r  c  d] 


(4.6.8) . 


The  individual  elements  in  F  can  be  determined  by  performing  the 

dU 

integration  from  0  to  w  of  the  expression  for  as  presented  in 


section  4.4.1.  Similarly  ,  the  individual  elements  in  D  can  be 

3W 

determined  by  performing  the  integration  of  the  expression  for  as 


presented  in  section  4.4.2.  To  perform  the  above  integrations,  as  well 
as  inverting  the  F  matrix  and  solving  for  the  x  's,  a  program  was 
written.  This  "Toroid  Program"  is  listed  in  Appendix  B. 
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The  "Toroid  Program"  was  run  to  compare  the  results  with  those 
obtained  vising  B0S0R4.  As  can  be  seen  from  Figures  4.4  and  4.5,  the 
results  obtained  for  the  displacements  are  comparable  to  those  obtained 
from  B0S0R4  and  are  of  the  same  general  shape.  Additionally,  when  the 
internal  energy  of  the  results  are  compared,  the  "Toroid  Program"  is 
only  4.3%  greater  than  that  of  the  B0S0R4  solution.  This  difference  can 
be  accounted  for  by  the  selection  of  the  functions  chosen  to  represent  u 
and  v.  The  functions  chosen  for  u  and  v  only  model  the  displacements 
and  do  not  model  the  change  in  displacements  or  the  curvature  of  the 
toroid. 

4.7  Forces  and  Moments: 

The  next  step  was  to  look  at  the  predictions  for  the  forces,  and 
,  and  the  moments,  and  ,  and  compare  them  with  the  results  from 
B0S0R4 . 

The  forces  and  moments  predicted  by  the  "Toroid  Program"  did  not 
corollate  well  with  the  results  obtained  from  B0S0R4.  As  can  be  seen 
from  Figures  4.6  and  4.7,  the  average  of  the  predicted  forces  is 
approximately  the  same  as  the  forces  obtained  from  B0S0R4.  However,  as 
can  be  seen  from  Figures  4.8  through  4.11,  the  predicted  moments  were 
several  orders  of  magnitude  less  than  those  obtained  from  B0S0R4. 


i 
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NORMAL  DISPLACEMENT 


Figure  4.4 

Normel  Displacement  (u) 


Figure  4 . 5 

Tangential  Displacement  (v) 
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Figure  4.6 
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Figure  4.7 
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Figure  4.8 
U,  vs  6  (B0S0R4) 
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Figure  4 . 10 
M,  v*  8  (B0S0R4) 
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Figure  4.11 

M,  v«  6  (Toroid  Program) 
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CHAPTER  5 


CONCLUSIONS 

As  was  shown  in  Chapter  4,  using  the  assumed  functions  for  the 
displacements  in  the  energy  method  produced  results  for  the 
displacements  that  were  consistent  with  those  obtained  from  B0S0R4. 
However,  it  was  also  shown  that  the  "Toroid  Program"  was  very  poor  at 
predicting  the  forces  and  moments  of  the  toroid. 

The  differences  can  be  explained  by  looking  at  what  the  assumed 
functions  used  for  the  displacements  represents.  The  assumed  functions 
for  u  and  v  are  fairly  good  models  of  the  displacements  but  were  not 
modeled  to  meet  the  slope  or  curvature  of  the  displacements  except  at 
the  end  points,  where  8  -  0  and  8  -  ir.  This  resulted  in  good 
predictions  for  u  and  v  but  inaccurate  predictions  for  u'  ,  u’  '  and  v'  , 
all  of  which  are  contained  in  the  expressions  for  the  forces  and 
moments . 

Some  thought  has  been  given  to  what  other  assumed  functions  might 
be  used  to  model  the  displacements  that  would  more  accurately  predict 
the  forces  and  moments.  It  is  believed  that  a  complete  Fourier  series 
representation,  instead  of  just  an  odd  or  even  series,  might  better 
predict  the  slope  and  curvature  of  the  displacements  and  therefore  the 
forces  and  moments.  To  ensure  that  the  complete  Fourier  series  meets 
the  boundary  conditions,  Lagrange  multipliers  would  be  required  to 
describe  the  constraint  energy  (See  section  4.2.3).  This  is  an  area 
that  is  left  for  further  investigation.  If  a  complete  Fourier  Series 


representation  does  not  adequately  predict  the  forces  and  moments  the 
next  step  would  be  to  write  some  kind  of  finite  element  or  finite 
difference  program  to  solve  this  problem. 

Since  the  program  written  does  model  the  displacements  and  produces 
results  that  are  consistent  with  those  obtained  from  B0S0R4,  the  "Toroid 
Program"  can  be  used  as  a  preliminary  design  tool  without  having  to 
resort  to  a  very  complex  and  expensive  computer  code  such  as  B0S0R4 . 
Unfortunately  the  "Toroid  Program"  is  restricted  to  analyzing  only 
toroids  under  hydrostatic  loading.  Complexities  in  describing  the 
geometry  of  a  general  shell  of  revolution  will  quickly  lead  the  designer 
to  the  more  complex  computer  codes  like  BOSOR4. 

It  is  also  believed  that  the  program  written  could  be  modified  to 
incorporate  nonlinear  terms  and  then  be  used  to  predict  buckling  loads. 
This  is  again  an  area  that  will  be  left  for  further  investigation. 


* 
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APPENDIX  A. 
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CURVE  FITTING  PROGRAMS 


The  programs  listed  were  used  in  conjunction  with  the  out  put  from 
B0SQR4  to  determine  the  coefficients  of  the  assumed  Fourier  Series. 
Both  programs  are  written  in  Fortran  77. 

Program  1  determines  ten  coefficients  for  an  odd  Fourier  Series 
that  is  symmetric  about  8  equal  to  to  fit  the  data  from  B0S0R4  for 
the  normal  displacement  of  the  shell  meridian. 

Program  2  determines  nine  coefficients  for  an  even  Fourier  Series 
that  is  anti -symmetric  about  8  equal  to  n,  to  fit  the  data  from  B0S0R4 
for  the  tangential  displacement  of  the  shell  meridian. 


PROGRAM  1  ULS.HPC1 

C  CURVE  FITTING  PROGRAM 
C  LEAST  SQUARES  CURVE  FITTING 
C234567 

INTEGER  I.J.M.N 

REAL  X(100) ,Y(100) , F( 100 , 10) , FT(10 , 100) , A(10 , 11) . B(10) 
REAL  FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 , F10 , C ( 10) 

EXTERNAL  FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 , F10 
N  -  100 
M  -  10 

C  READ  X-Y  VALUES  OF  DATA  POINTS 
DO  10  I-1,N 


READ  (UNIT  -  2 , FMT  -  *)  X(I) 
READ  (UNIT  -  3, FMT  -  *)  Y(I) 


10  CONTINUE 
C  GENERATE  THE  F  MATRIX 


20 

C 


I 

30 

C 

C 

I 

c 

c 


40 

C 

C 


C 


DO  20  I-l.N 

F(I.l)  -  F1(X(I)) 

F(1 , 2)  -  F2(X(I) ) 

F(I , 3)  -  F3(X(I) ) 

F(I,4)  -  F4(X(I) ) 

F(I , 5)  -  F5(X(I) ) 

F(I , 6)  -  F6 (X(I) ) 

F(I , 7)  -  F7 (X(I) ) 

F(I , 8)  -  F8(X(I) ) 

F(I , 9)  -  F9(X(I) ) 

F(I , 10)  -  F10(X(I) ) 

CONTINUE 

GENERATE  THE  TRANSPOSE  OF  THE  F  MATRIX 
DO  30  I-l.N 
DO  30  J-l.M 

FT( J , I )  -  F(I , J) 

CONTINUE 

DETERMINE  COEFFICIENT  MATRIX  A  OF  SIMULTANEOUS 
EQUATION  SYSTEM 

CALL  MATMPY (FT.F.A.M.N.M) 

DETERMINE  COLUMN  OF  CONSTANTS  FOR  SIMULTANEOUS 
EQUATION  SYSTEM 

CALL  MATMPY (FT,Y,B,M,N,1) 

DO  40  I-l.M 
A ( I , M+l )  -  B(I) 

CONTINUE 

DETERMINE  A(n)  VALUES  BY  SOLVING  SIMULTANEOUS  EQUATIONS 
USING  CHOLEiKY  METHOD 
MP1  -M+l 
CALL  CHLSKY(A , M, MP1 , C) 

WRITE  OUT  THE  A(n)  VALUES 
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WRITE  (UNIT  -  4 , FMT  -  *) (I . C(I) , 1-1 ,M) 

END 

C 

C  DETERMINES  MATRIX  C  AS  PRODUCT  OF  A  AND  B  MATRICES 
SUBROUTINE  MATMPY (A , B , C , M , N , L) 

REAL  A(M,N) ,B(N,L) ,C(M,L) 

DO  10  I-l.M 
DO  10  J-l.L 
C(I,J)  -  0. 

DO  10  K-l.N 

C(I,J)  -  C(I,J)+A(I,K)*B(K,J) 

10  CONTINUE 
END 

C  DEFINE  THE  FUNCTIONS 
C 

C  DEFINE  THE  A(0)  FUNCTION 
REAL  FUNCTION  F1(X) 

REAL  X 
FI  -  1.0 
END 

C  DEFINE  THE  A(l)  FUNCTION 
REAL  FUNCTION  F2(X) 

REAL  X 
F2  -  COS(X) 

END 

C  DEFINE  THE  A(2)  FUNCTION 
REAL  FUNCTION  F3(X) 

REAL  X 

F3  -  COS (2 .  *X) 

END 

C  DEFINE  THE  A(3)  FUNCTION 
REAL  FUNCTION  F4(X) 

REAL  X 
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F4  -  COS (3 .  *X) 

END 

C  DEFINE  THE  A(4)  FUNCTION 
REAL  FUNCTION  F5(X) 
REAL  X 

F5  -  COS (4 . *X) 

END 

C  DEFINE  THE  A(5)  FUNCTION 
REAL  FUNCTION  F6(X) 
REAL  X 

F6  -  COS (5 . *X) 

END 

C  DEFINE  THE  A(6)  FUNCTION 
REAL  FUNCTION  F7(X) 
REAL  X 

F7  -  COS (6 . *X) 

END 

C  DEFINE  THE  A(7)  FUNCTION 
REAL  FUNCTION  F8(X) 
REAL  X 

F8  -  COS ( 7 . *X) 

END 

C  DEFINE  THE  A(8)  FUNCTION 
REAL  FUNCTION  F9(X) 
REAL  X 

F9  -  COS (8 . *X) 

END 

C  DEFINE  THE  A(9)  FUNCTION 
REAL  FUNCTION  F10(X) 
REAL  X 

F10  -  COS ( 9 . *X) 

END 

C 
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SUBROUTINE  CHLSKY(A,N,M,X) 

REAL  A(10 , 11) ,X(10) 

C  CALCULATE  FIRST  ROW  OF  UPPER  UNIT  TRIANGULAR  MATRIX 
DO  10  J-2.M 
A(  1 , J )  -  A(1,J)/A(1,1) 

10  CONTINUE 

C  CALCULATE  OTHER  ELEMENTS  OF  U  AND  L  MATRICES 
DO  60  1-2, N 
J  -  I 

DO  30  II-J.N 
SUM  -  0. 

JM1  -  J-l 

DO  20  K-l.JMl 

SUM  -  SUM+A(II,K)*A(K,J) 

20  CONTINUE 

A(  II ,  J )  -  A( II , J ) -  SUM 
30  CONTINUE 
IP1  -  1+1 
DO  50  JJ-IPl.M 
SUM  -  0, 

IM1  -  1-1 

DO  40  K-l.IMI 

SUM  -  SUM+A(I ,K)*A(K,JJ) 

40  CONTINUE 

A(I.JJ)  -  (A(I , JJ) -SUM) /A (I , I) 

.50  CONTINUE 
60  CONTINUE 

C  SOLVE  FOR  X(I)  BY  BACK  SUBSTITUTION 
X(N)  -  A(N , N+l) 

L  -  N-l 
DO  80  NN-l.L 
SUM  -  0. 

I  -  N-NN 
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IP1  -  1+1 
DO  70  J-IPl.N 
SUM  -  SUM+A (I,J)*X(J) 
70  CONTINUE 

X(I)  -  A(I ,M) -SUM 
80  CONTINUE 
END 


PROGRAM  2  LISTING: 


C  CURVE  FITTING 

C  LEAST  SQUARES  CURVE  FITTING 

C234567 

INTEGER  I , J ,M,N 

REAL  X(100) , Y(IOO) , F( 100 ,9) , FT( 9 , 100) , A(9 , 10) , B( 9) 
REAL  FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 , C(9) 

EXTERNAL  FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 
N  -  100 
M  -  9 

C  READ  X-Y  VALUES  OF  DATA  POINTS 
DO  10  I-l.N 

READ  (UNIT  -  2 , FMT  -  *)  X(I) 

READ  (UNIT  -  5, FMT  -  *)  Y(I) 

10  CONTINUE 
C  GENERATE  THE  F  MATRIX 
DO  20  I-l.N 

F(I , 1)  -  F1(X(I) ) 

F(I , 2)  -  F2(X(I)) 

F(I , 3)  -  F3(X(I) ) 

F(I ,4)  -  F4(X(I) ) 

F(1 , 5)  -  F5(X(I) ) 

F(I , 6)  -  F6(X(I) ) 

F(1 , 7)  -  F7 (X(I) ) 

F(I , 8)  -  F8(X(I) ) 

F(I , 9)  -  F9(X(I) ) 

20  CONTINUE 

C  GENERATE  THE  TRANSPOSE  OF  THE  F  MATRIX 
DO  30  I-l.N 
DO  30  J-l.M 

FT(J , I)  -  F( I , J) 

30  CONTINUE 
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C  DETERMINE  COEFFICIENT  MATRIX  A  OF  SIMULTANEOUS 
C  EQUATION  SYSTEM 

CALL  MATMPY (FT.F.A.M.N.M) 

C  DETERMINE  COLUMN  OF  CONSTANTS  FOR  SIMULTANEOUS 
C  EQUATION  SYSTEM 

CALL  MATMPY (FT,Y,B,M,N,1) 

DO  40  I-1,M 
A(I ,M+1)  -  B ( I ) 

40  CONTINUE 

C  DETERMINE  B(n)  VALUES  BY  SOLVING  SIMULTANEOUS  EQUATIONS 
C  USING  CHOLESKY  METHOD 
MP1  -  M  +  1 
CALL  CHLSKY(A,M,MP1,C) 

C  WRITE  OUT  THE  B(n)  VALUES 

WRITE  (UNIT  -  4 , FMT  -  *) (I , C(I) , 1-1 ,M) 

END 

C 

C  DETERMINES  MATRIX  C  AS  PRODUCT  OF  A  AND  B  MATRICES 
SUBROUTINE  MATMPY( A , B , C , M , N , L) 

REAL  A(M,N) ,B(N,L) ,C(M,L) 

DO  10  I— 1 , M 
DO  10  J-l , L 
C(I,J)  -  0. 

DO  10  K-1,N 

C(I,J)  -  C(IIJ)+A(I,K)*B(K,J) 

10  CONTINUE 
END 
C 

C  DEFINE  THE  FUNCTIONS 

C 

C  DEFINE  THE  B(l)  FUNCTION 
REAL  FUNCTION  F1(X) 

REAL  X 
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FI  -  SIN(X) 

END 

C  DEFINE  THE  B(2)  FUNCTION 
REAL  FUNCTION  F2(X) 
REAL  X 

F2  -  SIN(2.*X) 

END 

C  DEFINE  THE  B(3)  FUNCTION 
REAL  FUNCTION  F3(X) 
REAL  X 

F3  -  SIN(3.*X) 

END 

C  DEFINE  THE  B(4)  FUNCTION 
REAL  FUNCTION  F4(X) 
REAL  X 

F4  -  SIN(4.*X) 

END 

C  DEFINE  THE  B(5)  FUNCTION 
REAL  FUNCTION  F5(X) 
REAL  X 

F5  -  SIN(5.*X) 

END 

C  DEFINE  THE  B(6)  FUNCTION 
REAL  FUNCTION  F6(X) 
REAL  X 

F6  -  SIN(6.*X) 

END 

C  DEFINE  THE  B(7)  FUNCTION 
REAL  FUNCTION  F7(X) 
REAL  X 

F7  -  SIN(7.*X) 

END 

C  DEFINE  THE  B(8)  FUNCTION 


58 


REAL  FUNCTION  F8(X) 

REAL  X 

F8  -  SIN(8 . *X) 

END 

C  DEFINE  THE  B(9)  FUNCTION 
REAL  FUNCTION  F9(X) 

REAL  X 

F9  -  SIN(9.*X) 

END 

C 

SUBROUTINE  CHLSKY(A,N ,M,X) 

REAL  A(9 , 10) ,X(9) 

C  CALCULATE  FIRST  ROW  OF  UPPER  UNIT  TRIANGULAR  MATRIX 
DO  10  J-2 ,M 
A(  1 , J )  -  A(1 , J)/A(l , 1) 

10  CONTINUE 

C  CALCULATE  OTHER  ELEMENTS  OF  U  AND  L  MATRICES 
DO  60  1-2, N 
J  -  I 

DO  30  II-J.N 
SUM  -  0. 

JM1  -  J-l 

DO  20  K-l.JMl 

SUM  -  SUM+A(II ,K)*A(K, J) 

20  CONTINUE 

A( II ,  J )  -  A(II , J) -SUM 
30  CONTINUE 
IP1  -  1+1 
DO  50  JJ-IPl.M 
SUM  -  0. 

IM1  -  1-1 

DO  40  K-l.IMl 

SUM  -  SUM+A(I,K)*A(K,JJ) 


40  CONTINUE 

A(I,JJ)  -  (A(I , JJ) -SUM)/A(I , I) 

50  CONTINUE 
60  CONTINUE 

C  SOLVE  FOR  X(I)  BY  BACK  SUBSTITUTION 
X(N)  -  A(N ,N+1) 

L  -  N-l 
DO  80  NN-l.L 
SUM  -  0. 

I  -  N-NN 

IPX  -  1+1 

DO  70  J-IPl.N 

SUM  -  SUM+A(I , J)*X(J) 

70  CONTINUE 

X(I)  -  A(I ,M) -SUM 
80  CONTINUE 
END 
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APPENDIX  B. 

TOROID  PROGRAM 


This  program  is  set  up  to  determine  the  ten  coefficients  of  the 
even  Fourier  Series  for  the  normal  displacement,  (u)  ,  of  the  toroid  and 
the  nine  coefficients  of  the  odd  Fourier  Series  for  the  tangential 
displacement,  (v) ,  of  the  toroid. 

To  determine  the  shape  and  magnitude  of  the  normal  and  tangential 
displacement  of  the  toroid  the  following  equations  are  used: 

I 

9 

u  -  7  a  cos(nS) 

n 

n  •  0 

9 

v  -  Y  b  sin(n0) 
n 

n  » 1 

The  program  was  written  in  Fortran  77 .  The  subroutines  used  were 
adopted  from  those  presented  by  Dyck,  Lawson  and  Smith  (1984) . 

PROGRAM  LISTING: 

C  TOROID  PROGRAM 
INTEGER  I ,N 

REAL  F(19 , 19) ,D(19) ,C(19) 

REAL  E , V , RB , RL , H , P , Al , A2 , A3 , A 
LOGICAL  ERROR 

EXTERNAL  FA, FAl , FB , FBI , SIMP 

PRINT  *,'This  Program  will  calculate  the  coefficients' 

PRINT  *,'of  the  Fourier  Series  for  the  Normal  and' 
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PRINT  *, 'Tangential  displacements  of  a  thin,  circular' 

PRINT  toroidal  shell  when  loaded  by  uniform' 

PRINT  *, 'external  hydrostatic  pressure.’ 

PRINT  * , ’  ' 

PRINT  *, 'Enter  Youngs  Modulus,  E(psi)  -' 

READ  * , E 

PRINT  *, 'Enter  Poissons  Ratio,  V  -' 

READ  *,V 

PRINT  *, 'Enter  Toroid  radius,  R(in)  -' 

READ  * , RB 

PRINT  *, 'Enter  Circular  Section  radius,  r(in) 

READ  * ,RL 

PRINT  *, 'Enter  Shell  thickness,  h(in)  -' 

READ  * ,H 

PRINT  *, 'Enter  hydrostatic  pressure,  P(psi)  -' 

READ  *,P 

PI  -  4 . 0*ATAN (1.0) 

N  -  19 
A  -  RB/RL 

Al-2.0*PI*E*H/(1. -V**2) 

A2--2.0*PI*P*RL**2 

A3-( -PI*E*(H**3) )/(6 ,0*(RL**2)*(1 .0-V**2) ) 

C  GENERATE  THE  D  CULUMN 
D(1)-A2*PI*A 
D(2)-A2*PI/2.0 
DO  30  1-3, N 
D(I)-0 .0 
30  CONTINUE 
C  GENERATE  THE  F  MATRIX 
DO  40  I-l.N 

PRINT  *,' INTEGRATING  MATRIX  COLUMN', I 

F(1 , I)  -  Al*SIMP(FA , PI, V, A, 1,0. 0)+A3*SIMP(FAl , PI, V, A, 1,0.0) 
F(2 , 1)  -  Al*SIMP(FA , PI, V, A, 1,1. 0)+A3*SIMP( FAl , PI, V, A, 1, 1.0) 
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F(3 , 1)  -  A1*SIMP(FA, PI ,V, A, 1 , 2 . 0)+A3*SIMP(FAl , PI, V, A, 1,2.0) 
F(4, I)  -  Al*SIMP(FA, PI , V, A, I , 3 . 0)+A3*SIMP(FAl ,PI,V,A, 1,3.0) 
F(5 , I)  -  Al*SIMP(FA, PI , V, A, I ,4 ,0)+A3*SIMP(FAl ,PI,VfA, 1,4.0) 
F(6,I)  -  Al*SIMP(FA, PI , V, A, I , 5 . 0)+A3*SIMP(FAl , PI, V, A, I, 5.0) 
F(7 , I)  -  Al*SIMP(FA, PI ,V, A, I , 6 . 0)+A3*SIHP(FAl , PI, V, A, I, 6.0) 
F(8 , I)  -  Al*SIMP(FA, PI , V,A, I , 7 . 0)+A3*SIMP(FAl , PI, V, A, 1,7.0) 
F(9 , I)  -  Al*SIMP(FA, PI , V, A, I , 8 . 0)+A3*SIMP(FAl , PI, V, A, I, 8.0) 
F(I0, I)  -  Al*SIMP(FA,PI ,V, A, I ,9 .0)+A3*SIMP(FAl , PI,V,A,I,9.0) 
F(ll , I)  -  A1*SIMP(FB,PI,V,A,I,1.0)+A3*SIMP(FB1,PI,V,A,I11.0) 
F(12 , I)  -  Al*SIMP(FB , PI ,V, A, I , 2 . 0)+A3*SIMP(FBl ,PI,V,A,I,2.0) 
F(13 , I)  -  Al*SIMP(FB , PI ,V, A, I , 3 . 0)+A3*SIMP(FBl , PI, V, A, I, 3.0) 
F(14, I)  -  Al*SIMP(FB , PI , V, A, I , 4 . 0)+A3*SIMP(FBl , PI, V, A, I, 4.0) 
F(15 , I)  -  Al*SIMP(FB , PI , V , A, I , 5 . 0)+A3*SIMP(FBl , PI, V, A, 1,5.0) 
F(16 , 1)  -  Al*SIMP(FB , PI , V , A, 1 , 6 . 0)+A3*SIMP(FBl , PI, V, A, 1,6.0) 
F(17 , I)  -  A1*SIMP(FB , PI ,V, A, I , 7 . 0)+A3*SIMP(FBl , PI, V, A, I, 7.0) 
F(18 , I)  -  Al*SIMP(FB,PI,V,A,I,8.0)+A3*SIMP(FBl,PI,V>A,I,8.0) 
F(19 , 1)  -  Al*SIMP(FB , PI , V, A, 1 , 9 . 0)+A3*SIMP(FBl ,PI,V,A,I,9.0) 
^0  CONTINUE 

C  CONVERT  THE  LINEAR  SYSTEM  TO  A  TRIANGULAR  SYSTEM 
CALL  GAUSS (F,D,N, ERROR) 

IF  (ERROR)  THEN 

PRINT  *, 'MATRIX  GENERATED  IS  SINGULAR,' 

PRINT  *, 'SOLUTION  IS  NOT  POSSIBLE.' 


C  SOLVE  THE  TRIANGULAR  LINEAR  SYSTEM 
CALL  BSOLVE(F,D,N,C) 

PRINT  * , ’ THE  FOLLOWING  ARE  THE  10  COEFFICIENTS  FOR  THE' 


PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 


'NORMAL  DISPLACEMENT  OF  THE  TOROID.’ 


'USE  IN  THE  SERIES;  An*COS (n*theta/r) ' 
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PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 


(I-l.C(I), 1-1,10) 

#  I 

'THE  FOLLOWING  ARE  THE  9  COEFFICIENTS  FOR  THE’ 
'TANGENTIAL  DISPLACEMENT  OF  THE  TOROID.’ 

t  $ 

'USE  IN  THE  SERIES;  Bn*SIN(n*theta/r) ' 


Bn 


+ 


Bn' 


PRINT  *,(1-10, C(I), 1-11,19) 

END  IF 

END 


C 

C  THIS  SUBPROGRAM  FUNCTION  INTEGRATES  THE  MATRIX  FUNCTIONS 
C  DEFINED  BY  USING  SIMPSON'S  RULE  AS  AN  APPROXIMATION 
REAL  FUNCTION  SIMP(F,PI ,V,A,N,Z) 

REAL  PI , F , H , SUMEVN , SUMODD , X , B 

INTEGER  I 

EXTERNAL  F 

B-0.0 

H-PI/100. 

SUMEVN-0 . 0 
SUMODD— F (A , V ,N , Z , H) 

DO  100  1-1,49 
X-2 . *I*H 

SUMEVN-SUMEVN+F( A , V , N , Z , X) 

SUMODD— SUMODD+F( A , V , N , Z , X+H) 

100  CONTINUE 


SIMP-(H/3 . )*(F(A,V,N,Z,B)+4. *SUMODD+2 . *SUMEVN+F (A,V,N,Z,PI)) 
RETURN 


C  DEFINE  THE  MATRIX  FUNCTIONS 
C  CALCULATE  ROWS  1  THROUGH  10 
C  FUNCTION  FA  REPRESENTS  MEMBRANE  ENERGY 
REAL  FUNCTION  FA(A,V, I ,Z,X) 

REAL  Q 

IF  (I  .LT.  11)  THEN 

Q-I-l 

FA-( A+COS (X) ) *C0S (Q*X)*COS ( Z*X) +2 . 0*V*C0S (X) *COS (Q*X) * 
+  COS(Z*X)+(COS(X)**2)*COS (Q*X)*COS(Z*X)/( A+COS (X) ) 

ELSE 

Q—I-10 

FA— (A+COS (X) )*Q*COS (Q*X)*C0S ( Z*X)+V*(C0S (X)*Q*COS (Q 
+  *X)*COS(Z*X) -SIN(X)*SIN(Q*X)*COS(Z*X) ) -COS (X)*SI 

+  N (X) *SIN (Q*X) *COS ( Z*X) / (A+COS (X) ) 

ENDIF 

END 

C  FUNCTION  FAl  REPRESENTS  BENDING  ENERGY 
REAL  FUNCTION  FA1(A,V, I , Z,X) 

REAL  Q 

IF  (I  .LT.  11)  THEN 

Q-l-1 

FAl- (A+COS (X) ) * ( (Q*Z)**2 ) *COS (Q*X) *COS ( Z*X) -  V*S IN (X 
+  )*(Q*(Z**2)*SIN(Q*X)*COS(Z*X)+(Q**2)*Z*COS(Q*X)* 

+  SIN ( Z*X) ) + ( SIN (X) **2 ) *Q*Z*SIN ( Q*X) *SIN ( Z*X) / ( A+C 

+  OS(X)) 

ELSE 

Q— I - 10 

FAl— (A+COS (X) )*Q*(Z**2)*COS (Q*X)*COS(Z*X) -V*SIN(X)* 

+  ( Q*Z*COS ( Q*X ) *S IN ( Z*X )+(Z**2)*SIN( Q*X ) *COS ( Z*X ) ) 

+  +( S IN (X)**2)*Z*SIN(Q*X)*SIN(Z*X)/( A+COS (X) ) 


C  CALCULATE  ROWS  11  THROUGH  19 
C  FUNCTION  FB  REPRESENTS  MEMBRANE  ENERGY 
REAL  FUNCTION  FB(A, V, I , Z ,X) 

REAL  Q 

IF  (I  .LT.  11)  THEN 
Q-I-l 

FB- ( A+COS (X) ) *Z*COS (Q*X) *COS ( Z*X) +V* ( COS (X) *Z*COS (Q 
+  *X)*COS(Z*X) -SIN(X)*COS(Q*X)*SIN(Z*X) ) -COS (X)*SI 

+  N(X)*COS(Q*X) *SIN(Z*X)/ (A+COS (X) ) 

ELSE 

Q— 1-10 

FB-( A+COS (X) ) *Q*Z*COS (Q*X) *COS ( Z*X) *  V*SIN  (X) * ( Z*SIN 
+  (Q*X) *COS ( Z*X) +Q+COS ( Q*X) *S IN ( Z*X) ) + ( S IN (X) **2 ) * 

+  SIN(Q*X)*SIN(Z*X)/(A+COS(X) ) 

END  IF 
END 

C  FUNCTION  FBI  REPRESENTS  BENDING  ENERGY 
REAL  FUNCTION  FBI (A , V , I , Z ,X) 

REAL  Q 

IF  (I  .LT.  11)  THEN 
Q-I-l 

FBI- (A+COS (X) ) *Z* (Q**2) *COS (Q*X) *COS ( Z*X) - V*S IN (X) * 

+  (Q*Z*SIN(Q*X)*COS (Z*X)+(Q**2)*COS (Q*X)*SIN(Z*X) ) 

+  +( SIN (X)**2)*Q*S IN (Q*X)*SIN(Z*X)/( A+COS (X) ) 

ELSE 

Q-I-10 

FBI- (A+COS (X) ) *Z*Q*COS (Q*X) *COS ( Z*X) -V*SIN(X)*( Z*S I 
+  N(Q*X)*COS ( Z*X) +Q*COS (Q*X) *SIN ( Z*X) ) + ( S IN ( X) **2 ) 

+  *SIN(Q*X)*SIN(Z*X)/(A+COS (X) ) 

END  IF 
END 
C 

C  GAUSSIAN  ELIMINATION  MODULE 
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SUBROUTINE  GAUSS ( A, B,M, ERROR) 
INTEGER  M, I , J ,K, INDEX 
REAL  A(M,M) ,B(M) .PIVOT, TEMP, RATIO 
LOGICAL  ERROR 
ERROR- . FALSE . 

DO  100  I-l.M 
INDEX-I 

PIVOT— ABS (A( I , I ) ) 

DO  200  J-I+l.M 

IF  (ABS (A( J , I) )  .GT.  PIVOT)  THEN 
PIVOT— ABS (A(J , I) ) 

INDEX-J 
ENDIF 
200  CONTINUE 

IF  (INDEX  .GT.  I)  THEN 
DO  400  K-I.M 
TEMP— A( I , K) 

A(I ,K)— A( INDEX, K) 

A( INDEX, K)-TEMP 
400  CONTINUE 
TEMP-B ( I ) 

B(I)-B(INDEX) 

B( INDEX) -TEMP 
ENDIF 

IF  (PIVOT  . EQ.  0.0)  THEN 
ERROR  -  .TRUE. 

ELSE 

DO  300  J-I+I.M 

RATIO— A(J,I)/A(I,I) 

DO  500  K-I+l.M 

A(J,K)-A(J,K)-A(I,K)*RATIO 


500 


CONTINUE 

B(J)-B(J) -B(I)*RATI0 


300  CONTINUE 
END  IF 

100  CONTINUE 
END 
C 

C  SUBROUTINE  TO  SOLVE  TRIANGULAR  LINEAR  SYSTEM 
SUBROUTINE  BSOLVE(A, B ,M, Z) 

INTEGER  M , I , J 
REAL  A(M,M) ,B(M) ,Z(M) , SUM 
DO  200  I-M.l, -1 
SUM-B(I) 

DO  100  J-I+l.M 

SUM— SUM -A(I,J)*Z(J) 

100  CONTINUE 

Z ( I ) —SUM/A (1,1) 

200  CONTINUE 


